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1.  INTRODUCTION 


This  report  describes  a  supplement  to  the  TIGER  computer  program. }  The 
latter  is  a  program  for  the  calculation  of  thermodynamic  properties  of 
nonideal  systems  of  specified  atomic  compositions  containing  gaseous,  liquid 
and  solid  phases  with  known  equations  of  state.  The  TIGER  program  version 
that  Is  documented  in  Ref.  1  has  options  to  use  for  the  gaseous  constituents 
either  the  ideal  equation  of  state  or  one  of  three  nonideal  equations.  The 
latter  three  are  called  BKW  (Becker,  Kistiakowsky  and  Wilson),  JC21  and  JCZ2 
(Jacobs,  Cowperthwaite  and  Zwisler) ,  respectively,  and  defined  in  Ref.  1.  The 
present  report  describes  a  modification  of  the  TIGER  program  which  permits  one 
to  also  use  a  fifth  gaseous  state  equation,  namely  ,v  a  truncated  virial 
equation  in  which  the  second  virial  coefficient  is  determined  from  lennard- 
Jones  (6-12)  potential  parameters*  The  added  equation  of  state  subroutine  has 
been  given  the  name  U612  in  order  to  indicate  its  relation  to  the  Lennard- 
Jones  potential.  The  principal  features  of  the  added  equation  of  state  are 
described  in  Section  2,  and  control  cards  that  activate  the  routine  LJ612  are 
discussed  in  Section  3.  Section  4  gives  a  sample  calculation.  Details  of  the 
derivation  of  the  pertinent  equations  and  formulas  are  given  in  Sections  5,  6 
and  7.  These  details  should  be  of  interest  for  users  of  the  U612  option  if 
the  basic  equations  are  modified,  for  instance,  for  numerical  expediency.  In 
order  to  facilitate  such  a  modification  the  subroutines  U612  and  UBS  (an 
auxiliary  routine)  are  listed  in  Appendices  A  and  B,  respectively,  with 
numerous  explanatory  comments.  These  comments  make  the  listing  of  U612  also 
useful  as  a  guide  for  the  preparation  of  other  nonideal  equations  of  state 
routines.  Section  8  contains  a  summary  and  conclusion. 

In  order  to  u«e  the  subroutine  U612  in  the  TIGER  program,  the  subroutine 
STATEG  rnst  be  changed  as  described  in  Ref.  1,  p.  III-C-34&.  That  change 
merely  increases  by  one  the  number  of  available  gaseous  equations  of  state  and 
includes  a  call  to  the  added  subroutine  U612  at  the  proper  place.  No  other 
changes  in  the  TIGER  program  are  needed. 


2.  PRINCIPAL  RESULTS 

The  modifications  of  the  TIGER  program  permit  one  to  use  the  following 
type  of  equation  of  state  for  gases: 


P  =  $(V,T) 

V 


(2.1) 


Ccwpertfiixiite  and  W.  H.  Zwisler,  "TIGER  Computer  Program  Documentation ,  r 
Stanford  Research  Institute  Publication  7206,  January  1573. 
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In  these  equations,  p  (Pa)  is  the  pressure,  R  *  8.3143  J/(K»mol)  is  the 
universal  gas  constant,  T  (k)  is  temperature,  V  (n^/raol)  is  the  molar 
volume,  v  is  the  imperfection  parameter,  B  (cr/mol)  is  the  second  virial 
coefficient,  and  C  (m^/aol^)  is  the  third  virial  coefficient.  The  formula 
(2.2)  for  $  is  called  a  truncated  virial  form  because  it  <"an  be  considered 
as  the  first#  three  terms  of  an  infinite  series  expansion  with  terms  of  the 
type  Aj»(T)/V^.  Ihe  second  virial  coefficient  B(T)  In  Eq.  (2.2)  is  computed 
assuming  that  the  intermolecular  forces  have  the  empirical  potential 
function 


$  (r) 


4  e  ((a/r)12  -  (o/r)6| 


(2.3) 


where  r  (m)  is  the  distance  between  molecules,  and  e  (J)  and  a  (a)  are 
parameters  characteristic  of  each  gas  (Ref.  2,  p.  32).  the  function  defined 
by  Lc .  (2.3)  is  called  the  Lennard-Jones  (6-12)  potential,  a  that  '»»’uc 
of  r  f,  r  which  $(r)  =  0,  and  e  is  the  maximum  energy  of  attraction  (or  der'ti 
of  th?  potential  well)  which  occurs  at  r  =  2*^0  =  1.12  o.  Typical  values 
of  <?  are  around  3.5*10~*°  ta.  The  parameter  e  enters  the  present 
calculations  only  in  the  ratio  e/k,  where  k  =  1 .38054* J/K  is  the 
Boltzmann  constant.  Typical  values  of  the  ratio  e  A .  are  around  300  K. 

For  a  single  gas  with  giver  o  and  e/k  the  second  virial  coefficient 
B(T)  is  calculated  as  follows.  vjfcrf.  2,  p.  162  ff./.  First,  one  defines  a 
dimensionless  reduced  temperature  T  v 


T  =  T  k/e 


(2.4) 


and  a  reduced  dimensionless  coefficient  B  by 


* 

B 


B/b  (c> 
o 


(2.5) 


where 


Zj.Q.  Hirschf elder,  C.P.  Curtiss  and  R.B.  Bird,  "Molecular  Theory  of  Cases 
and  Liquids,"  John  Wiley  and  Sons ,  Sen)  York,  1954. 


b  (o)  =  4  *  N  a3  •  (2.6) 

o  J 

V  JO  _»  * 

and  N  =  6. 02252*  10iJ  mol  Is  the  Avogadro  number.  B  is  for  a  given 
Jntermolecular  potential  function  of  Lennard-Jones  type  a  unique  function  of 
T.  It  can  be  calculated  by  a  series  expansion  and  it  is  tabulated  for  the 
Lennard-Jones  (6-12)  potential  in  Ref.  2,^p.  1114  ff.  Ihe  relation  between 
B,  T,  cf  and  e  is  in  terns  of  the  function  B  as  follows 


B(T)  =  b0(o).B(T  k/e)  .  (2.7) 

For  gas  mixtures  the  computation  of  B(T)  is  modified  such  that  in  Eq.  (2.7) 
one  uses  for  o  and  e  averages  of  the  parameters  of  pairs  of  individual 
constituents,  and  computes  the  final  B(T)  by  averaging  over  all  pairs. 
Details  of  the  calculation  of  B(.T)  are  given  in  Sections  5  and  6. 

A  corresponding  computation  of  the  third  virial  coefficient  C  for  gases 
with  Lennard-Jones  potential  is  described  on  Ref.  2,<  pp.  150  ff.,  170  ff., 
228  ff.  and  1119.  Ihe  algorithm  is  complicated  and  the  resulting  C  is  not 
accurate.  Therefore,  instead  of  using  the  Lennard-Jones  potential,  the 
third  virial  coefficient  computation  for  a  single  gas  is  for  the  present 
task  based  on  a  rigid  sphere  approximation  from  Ref.  2,  p.  157: 


C(o)  =  |  b2(o).0.Sl6  , 
o  o 


(2.8) 


where  bQ(cj)  is  given  by  Eq.  (2.6).  To  simplify  further,  the  third  virial 
coefficient  for  mixtures  of  gases  is  computed  by  a  simple  averaging  of  the 
individual  values  of  C,  described  in  Section  5.  Hirschfelder  et  al.^  p. 
153,  derive  a  more  complicated  fortaxla  but  considering  the  large 
uncertainties  of  the  result,  the  simpler  averaging  was  deemed  to  be 
adequate.  This  computation  of  the  third  virial  coefficient  was  suggested  by 
E.  Freedman  and  it  is  implemented  in  the  BLAKE  computer  program. ^  the 
factor  0.81^  is  empirical  and  suggested  in  Ref.  2,  p.  157.  It  may  be 
replaced  by  a  different  factor,  as  described  in  Section  3. 
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In  summary,  the  supplemented  TIGER  program  can  be  used  with  the 
equation  of  state  (2.1)  and  (2.2)  for  gases.  The  user  should  supply  for 
each  gaseous  constituent  of  the  system  the  proper  values  of  o  and  e  of  the 
Lennard-Jones  (6-12)  potential  parameters.  If  such  are  not  supplied  and 
also  not  already  stored  In  the  TIGER  file  of  material  constants,  then  the 
LJ612  routine  uses  the  following  default  values: 


a  =  3.5*10 


■10 


(2.9) 


e/k  =  300  K 


3.  USER’S  GUIDE  FOR  THE  SUBROUTINE  LJ612 

The  virial  equation  of  state  supplement  LJ612  of  the  TIGER  program  can 
be  activated  by  the  same  type  of  control  cards  as  the  other  four  equations 
of  state.  Particulars  of  input  data  and  order  are  provided  in  Ref.  1.  This 
section  supplements  that  reference,  particularly  its  Volume  IV,  ‘User’s 
Guide  of  the  TIGER  Computer  Program.” 

The  instruction  to  use  the  LJ612  equation  of  state  is  giv^n  by  a  GEOS- 
card  which  has  the  format 


GEOS,U612 


(see  Ref.  1,  IV-C-10-11).  The  card  instructs  the  TIGER  program  to  use  the 
U612  equation  of  state  for  all  computations  until  another  GEOS  card  is 
encountered  in  the  input.  Therefore,  it  should  be  placed  in  the  input  deck 
before  any  of  the  instruction  cards  that  initiate  a  computation,  that  is, 
before  POINT,  ISOLINE,  GRID,  EXPLOSION,  C-J  CONDITION  or  HUGONIOT  cards. 

The  empirical  factor  0.81  that  enters  Eq.  (2.8)  for  the  third  virial 
coefficient  can  be  changed  by  a  SET-card.  (See  Ref.  1,  IV-C-18).  The  card 
has  the  format 


SET,LJ612,SFACT,a 


where  “a"  is  the  number  that  replaces  the  value  0.81  in  Eq.  (2.8).  For 
instance,  if  a  =  0,  then  C  *  0  and  the  virial  Eq.  (2.2)  is  truncated  to  two 
terms.  The  SET-card  should  of  course  precede  all  the  calculation 
instruction  cards  listed  above.  Once  the  factor  has  been  set  then  it 
remains  in  effect  for  the  particular  TIGER  run  until  changed  by  another  SET- 


card. 


The  Lennard-Jones  (6-12)  potential  parameters  o  and  e/k  are  provided  by 
a  STG-card  for  each  gaseous  contituent.  These  cards  are  used  in  a  TIGER  run 
that  updates  (replaces)  the  TIGER  file  containing  material  constants.  The 
input  and  instruction  cards  for  such  a  run  are  desc-ibed  in  Ref.  I,  IV-C-24 
ff.  For  a  gaseous  constituent  with  the  designation  "name”  the  input 
consists  of  the  following  sequence  (see  Ref.  1,  IV-C-25): 


CONSTITUENT, name, GAS 
STR, name, GAS, 1 ,  •  •  • 
STR, name, GAS, 2,  .  .  . 
STR,name,GAS,3,  .  .  . 
STG,LJ612,name,o ,  e/k 


Other  STG-cards  specifying  constants  for  the  BKW,  JCZ2  or  JCZ3  equations  of 
state  may  follow  or  precede  the  STG,LJ612-card.  The  next  card  after  the 
STG-cards  is  another  CONSTITUENT-card.  The  value  of  a  must  be  entered  in 
10"10  n  (in  Angstroms)  and  the  value  of  e/k  should  be  expressed  in 
kelvins.  An  example  of  a  STG-card  is 


STG,LJ612,H20,2.79 ,542.5 


indicating  that  for  the  constituent  called  H20  the  Lennard-Jones  potential 
parameters  are  a  =  2.79* lO*1^  m  and  e/k  *  542.5  K.  If  either  o  or  e/k  is 
negative  or  zero  in  the  STG-card,  then  the  corresponding  default  value  (o  = 
3.5*  10“i0  n,  e/k  =  300  K)  is  stored  in  the  library  file  instead  of  the 
negative  or  zero  value  from  the  card.  If  the  library  file  does  not  contain 
values  of  the  Lennard-Jones  potential  parameters  for  a  constituent,  then 
default  values  are  generated  as  needed  by  the  subroutine  LJ612  at  run 
time.  (That  is,  if  the  LJ612  is  activated  by  a  GEOS ,LJ6 12-card  and  the 
constituent  is  part  of  the  system  that  is  evaluated. )  The  contents  of  the 
library  file  are  not  affected  by  such  computations. 

More  details  of  the  structure  and  operation  of  the  subroutine  U612  and 
the  auxiliary  routine  LJBS  are  provided  by  comments  in  the  listings  of  the 
routines  in  Appendices  A  and  B,  respectively. 


S 
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4.  EXAMPLE  OF  A  COMPUTATION 

We  chose  as  an  example  for  the  operation  of  the  LJ612  equation  of  state 
routine  the  same  problem  as  in  Ref.  1,  p.  IV-E-17  ff.  The  present  input 
differs  from  that  in  Ref.  1  only  by  the  GEOS-card.  The  output  is  given  in 
Appendix  C.  Comparison  of  the  output  with  that  of  Ref.  1  shows  two  types  of 
differences.  First,  one  observes  the  expected  difference  in  values  for 
temperature,  pressure  and  other  thermodynamic  variables,  because  a  different 
equation  of  state  was  used.  However,  in  addition  to  these  changes,  and,' 
possibly  masking  the  effect  of  the  equation  of  »tate,  there  is  also  a  change 
of  the  composition  of  the  gas.  This  change  is  not  a  consequence  of  the  new 
equation  of  state  but  rather  of  the  evolution  of  the  TIGER  library  data 
between  the  publication  of  Ref.  1  and  its  present  version. 


5.  DERIVATION  OF  PRINCIPAL  FORMULAS 


We  consider  an  equation  of  state  of  the  form 


P  *(V,T) 

V 


$  ( V ,  T )  =  1  rM+  <L  .  (5.2) 

V 

In  the  TIGER  program,  the  molar  volume  V  is  not  directly  used  as  a 
variable.  Instead,  it  may  be  calculated  from  the  following  relations  (Ref. 
1,  p.  I-B-16): 

,  M  M 

V=-5-v=--S-  =  -^-  (5.3) 

n  n  p  n  ^ 


where  n  is  the  number  of  moles,  V  is  the  gas  volume.  Mg  is  the  mass  of  the 
gas,  p  is  its  density  and  MQ  is  the  mass  of  the  mixture.  The  density  £  is 
defined  as  the  mass  of  the  mixture  divided  by  the  volume  of  the  gas.  The 
quantities  n,  MQ  aud  £  are  TIGER  variables  from  which  V  may  be  computed. 
Let  the  system  have  s  gaseous  constituents  with  mole  numbers  i=i,  ..., 
s.  Then 


n  “  l  n, 

1=1 


10 


The  ni  are  also  TIGER  variables.  The  equation  of  state  is  in  terms  of  these 
variables. 


p  =  R  T  n  r  J 


$  »  $(p,  T,  Mq,  n^,  n2,  ...»  ng) 

a  1  +--n  BCT,!^,  ...,ng)  +  (^-J2  n2  C(n;,  ...»  ng) 


The  virial  coefficients  B  and  C  depend  on  the  mole  numbers  n^  of  the 
individual  constituents  because  they  represent  the  non-ideal  effect  of  al-1  s 
constituents  and  the  nt  specify  the  quality  of  the  mixture.  Both 
coefficients  also  depend  on  the  individual  Lennard-Jones  potential 
parameters  and  e^,  In  a  form  that  will  be  shown  shortly. 

The  second  virial  coefficient  B  is  calculr^ed  for  a  system  of  gases  by 
the  formula  (Ref.  2,  p.  153) 


s  s 

;  — r  l  I  n  n.  B 
n2  1=1  1  =  1  '  J  ij 


In  this  formula  the  are  the  second  virial  coefficients  of  a  potential 
function  which  describes  the  interaction  between  molecules  of  constituent  i 
and  constituent  j.  They  are  functions  of  the  form 


Bij  =  Bi.(T,a1,o.,ei/k,£j/k)  , 


that  is,  they  depend  on  the  gas  temperature  T  and  on  the  Lennard-Jones 
potential  parameters  of  both  constituents.  The  calculation  of  the  B^  will 
be  described  in  Section  6. 

The  third  virial  coefficient  C  in  Eq.  (5.6)  is  computed  as  the  average 

C  =  -  y  n,  C.  (5.9) 

n  1=1  11 

of  the  third  virial  coeffients  of  the  individual  constituents.  The 
latter  are  assumed  to  be  independent  of  temperature  and  that  is. 


(5.18) 


5  In  $  _  J_  53 
5  n j  3  3 


+  Ic' 


+  £-i  V 


j-i 


Vu +  ‘k] 


n  C> 


Notice  that  the  first  two  expressions,  Eqs.  (5.16)  and  (5.17),^  are 
dimensionless*  whereas,  the  expression  (5.18)  has  the  dimension  1/nol. 
Another  dimensionless  quantity  is 


3  -  1 


,  n  53 
+  n  - — -j  7T 
on.J  p 


dp' 


(5.19) 


An  evaluation  of  the  Integral  (5.19)  yields  for  I*  the  formula 


7 1  “  2  M"  1  “jBIJ  +  I  D\  +  *C 

oj3l  J  J  O  o 


(5.20) 


Other  expressions  needed  by  the  TIGER  program  are  the  following  two 
dimensionless  derivatives  of  1*^: 


ffer  - *  w  * 2  r  l.  njBij +  f-§->2  "V  2 1&2*  Ic 

O  J-I  J  J  o  o 


and 


(5.21) 


ar.  ar.  ~  s  dB. 

t  „  T _ i  e  y  _  T _ L 

3  In  T  3T  H  / ,  V  dT 

o  j*l 


and  the  following  derivative  with  the  dimension  1/mol 

37r-2^-Bij  +  (M')2iIc  +  n(ci  +  cj)j  • 

j  o  J  o  J 


(5.22) 


(5.23) 


The  next  two  quantities  have  the  dimension  of  moles  and  are  defined  by 
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p  M 

.  _  r  _ O 

“  "  '  FT  p 


lo  ,  3p.  dp  9  T  3(ni)  d£ 

7 (p " T  ^  r  1  t_3 ~f 


/■,  M  T  ^ 

’  r  _a_  ijp  —  a(lt) 

eT"'-'  Ri?  ,_2  o'  “  3T 
o  3T 

An  evaluation  of  these  expressions  yields 


r  =  T  v 
*  M  -  B 


eT  “  M  (2  T-B  +  ^  -B  >  *  (5*27) 

*  o 

We  notice  that  the  quantity  e  that  is  defined  by  Eq.  (5.24)  or  (5.26)  is  an 
auxiliary  variable  for  the  TIGER  program  and  not  the  Lennard-Jones  potential 
parameter  e . 

In  order  to  evaluate  the  expressions  (5.15)  through  (5.18),  (5.20) 
through  (5.23),  (5.26)  and  (5.27),  one  needs  values  for  the  virial 
coefficients  Bjj  and  its  temperature  derivatives  which  enter  the  formulas  in 
the  combinations  T-dBjj/dT  and  T2»d2Bij/dT^,  and  the  virial  coefficients 
Cj.  The  computation  of  the  former  three  is  treated  in  the  next  section. 
The  computation  of  the  is  described  In  Section  7. 


6.  COMPUTATION  OP  SECOND  VIRIAL  COEFFICIENTS 


The  second  virial  coefficient  Bjj  in  Eq.  (5.7)  represents  the 
interaction  between  molecules  of  constituent  i  and  constituent  j.  It  may  be 
computed  as  the  virial  coefficient  of  a  pure  substance  with  an 
intermolecuiar  potential  that  is  an  empirical  combination  of  the  potential 
functions  of  both  constituents.  In  this  report,  we  use  the  following 
empirical  combining  rules  to  specify  the  constants  a  and  e . .  of  the 
mixture  potential  (Ref.  2,  p.  168,  Ref.  3,  p.  11):  ?  ^ 


°ij =  i  (ci  +<y 


-  I 

1 

4 


lii.fi.  1L13/2  '_l.  :i>'/2 

k  '‘rt  n  >  '‘k  k  > 


i 


We  notice  that  for  l*j  these  formulas  produce  0  *  <?i ,  and, 

consequently,  11  B^.  The  0^  and  Cj*  are  treated  as  the  Lennard-Jones 

(6-12)  occential  parameters  of  aJpure  substance  and  Bjj  is  computed  as  the 
second  virial  coefficient  of  that  substance.  This  is  done  by  defining  a 
dimensionless  reduced  temperature  J  and  confuting  a  corresponding 
dimensionless  virial  coefficient  B(T)  from  which  Bjj  then  can  be 

retrieved.  The  function  B(T)  is  uniquely  determined  by  the  constants  of  the 
Lennard-Jones  potential.  The  following  formulas  are  used  to  carry  out  this 
calculation: 

T  =  T  k/e  .  ,  (6.3) 


.  2-3 

b  =  ■=■  r.  No., 
oij  3  ij 


Bu  -  bolj  *B(« 


where  N  =  6. 02252* 102j  ool"1  Is  the  Avogadro  constant. 

the  function  B(T)  is  defined  by  the  following  infinite  series  (Ref.  2, 
p.  163): 


B(T)  -lb!  <-1/4  ‘  °/2> 


b  =  -  — r(-  i  +  I», 

D  4*  q!  4  L 


The  first  two  coefficients  in  the  series  (6.6)  are* * 


qF.  Jahrke  and  5.  ?ndet  "Tables  of  Higher  Functions B.G.  Teubner,  Leipzig 
2946. 


(6.8) 


bC  =  -  ^  r(-  -  n-  1.225  416  70 2  7 


and 


bi  -  r  r(4) 


/T 

2 


3.625  605  908 


(6.9) 


U1  other  coefficients  can  be  calculated  froa  b^  and  bj  by  the  recurrence 
formula 


b_.,  *  b  7^iT-7rT7T  ,  o  =  G,.  1,  2,  ... 
a+2  m  (m+1)  (m+2) 

We  rearrange  the  series  (6.6)  as  follows  for  aritncetical  expediency 


*A-  T  [b2m!^-“>  +  b2afl*T<-^-“)] 

a=G 


(6.10) 


(6.11) 


-  I  (A0'm+Al,o) 

a=0 

The  Aq.jjj  and  Aj  m  are  defined  by 

\>,o  ■  *•»  *  (-1/4)  • 

*1.0 -b, 


(6.12) 


(6.13) 


and  by  the  recurrence  formulas  for  a  =  0,1,2,... 


4a  -  1 


A0,c+1  =  *0.0  (2nd  1)  (2n+2)  *  * 


(6.14) 


4o  +  1 


1 


'M.nrH  Al,m  (2nH-2)  (2m*-3)  * 


(6.15) 


The  series  for  B(T)  converges  absolutely  for  all  |  t|  >  0.  Therefore,  one 
can  compute  the  derivatives  of  B(T)  by  termwise  differentiation  which 
pr  duces  the  following  formulas 


"M 


X 


Mb* 


T^=  I 


2*  oo 

t24|  -  I  [ 


(4mfl)  (4m4-5) 


(4nf3)  (4crf-7) 


Tne  evaluation  of  the  series  (6.11),  ^6.16)  and  (6.17)  is  done  in  the 
£ub£0>^tine  U£|* which  computes  for  given  T  the  corresponding  ^values  of 
B,  T  B*  ana  T  B"  .  The  series  converge  very  fast  T  is  large,  but 
require  the  confutation  of  a  large  number  o£  terms  if  T  is  small.  Typical 

TIGER  computations  3re  done  with  values  of  T  between  one  and  t^o^  In  that 
range  some  computation  time  can  be  saved  by  interpolating  B(T)  and  its 
derivatives  in  tables  instead  of  evaluating  the  series.  Such  tables  are 

published  in  Ref.  2,  p.  1144  ff..  Tables  I-B.  The  subroutine  UBS  uses 
the^e  tables  for  the  computation  of  B  and  its  derivatives  by  interpolation 
if  T  is  within  the  range  (0.75,  5.0).  If  T  *  5.0,  then  only  six  to  seven 
tercjs  of  the  series  have  to  be  calculated,  and  the  number  of  terns  decreases 
as  T  increases.  (The  end  criterion  for  the  series  calculation  in  UBS  is 

that  the  ^bsolute  value  of  the  last  term  is  less  than  10  . )  On  the  other 

hand,  if  T  <  0.75  then  the  table  interpolation  Is  not  very  accurate  due  to 
rapid  changes  of  the  functions.  Therefore,  within  the  range  (0.01,  0.75) 

again  the  series  formulas  are  used.  The  number  of  terms  that  have  to  be 

evaluated  is  350  at  T  =  0.01.*  However,  one  would  not  expect  typical  TIGER 
calculations  to  be  done  for  T  <  0.75.  For  T  <  0.01  an  error  stop  with  a 
corresponding  printed  message  is  programmed  into  U612. 

The  interpolation  in  the  tables  is  done  by  Heraite  interpolation 
formulas  which  make  use  of  the  availability  of  the  derivatives  of  the 
interpolated  functions.  Particulars  about  the  interpolation  formulas  are 
given  in  Appendix  D. 

In  a  limited  number  of  test  runs  the  computing  time  that  was  saved  by 
interpolation  was  found  to  be  of  the  order  of  two  percent. 


7.  COMPUTATION  OF  THIRD  VIRIAL  COEFFICIENTS 

The  third  virial  coefficient  C  is  computed  by  Eq.  (5.9)  as  an  average 
of  the  third  virial  coefficients  Cj  of  the  individual  constituents.  The 
latter  coefficients  we  compute  from  a  rigid  sphere  approximation.  According 
to  Ref.  2,  p.  157,  a  reasonable  value  for  in  case  of  high  temperature 
powder  gases  is 


irsww 


mol  *  is  the  Avogadro  constant  and  a  is  the  Leinard- 


where  K  =  6.02252* 10^ 

Jones  (6-12)  potential  parameter.  The  factor  0.81  in  Eq.  (7.1)  is  an 
approximation  to  the  theoretical  value  for  the  Lennard-Jones  potential  at 

60.  (|ee  Eqs.  (2.4)  and  (6.3)  for  the  definition  of  the  reduced 


^emperature  T. )  Theoretically  the  factor  nay  be  increased  to  about  0.99  for 
T  =  1.25  (as  can  be  seen  fron  Ref.  2,  pp.  171  and  1116  -  1117).  However, 
experimental  measurements  deviate  considerably  fron  the  theoretical  vaiae 
and,  therefore,  the  constant  0.81  likely  suffices  to  establish  a  correct 
order  of  magnitude  effect  of  the  third  virial  coefficients  on  TIGER 
calculations.  If  necessary,  the  factor  can  be  changed  for  any  particular 
calculations  by  using  the  SET-card.  (See  Section  3).  Particularly,  by 
setting  the  factor  equal  to  zero  and  repeating  the  calculation  one  can 
obtain  the  total  effect  of  the  third  virial  coefficients  on  the  calculated 
results. 


8.  SDKMARY  AND  CON CL OS I ON 


This  report  describes  two  new  subroutines  tnat  can  be  included  in  the 
TIGER  computer  program  and  enables  the  latter  to  do  computations  using  a 
truncated  virial  equation  of  state  for  the  gaseous  constituents.  The  second 
virial  coefficient  in  the  equation  is  based  on  Lennard-Jones  (6-12) 
intermolecular  potential,  and  the  third  virial  coefficient  is  based  on  a 
simplified  rigid  sphere  assumption.  The  new  equation  of  state  subroutine  is 
called  LJ£12.  It  uses  an  auxiliary  subroutine  LJBS  which  conputes  a 
function  B  and  its  derivatives  for  given  T.  The  latter  is  a  reduced 
dimensionless  temperature  and  B  is  a  dimensionless  second  virial  coefficient 
for  Lennard-Jones  (6-12)  potential  and  pure  substances.  The  U612 
subroutine  has  an  option  (SET-card)  to  multiply  the  third  virial  coefficient 
by  an  arbitrary  constant.  By  setting  the  constant  equal  to  zero  one  can 
obtain  an  estimate  of  the  significance  of  the  third  virial  coefficient  for 
any  particular  computation. 


The  TIGER  program  was  carefully  studied  during  the  programming  of  the 
additional  routines.  It  was  found  that  the  description  of  the  program  in 
Ref.  1  is  not  up  to  date  and,  particularly,  that  the  comments  in  the  BKW 
equation  of  state  routine  are  not  complete  and  sufficient  for  an  easy 
addition  of  new  equations  of  state.  Therefore,  tbe  new  subroutine  LJ612  was 
provided  with  numerous  comments  whicn  may  be  helpful  for  adding  another 
equation  of  state  to  the  TIGER  program.  We  also  notice  that  although  the 
TIGER  program  guide  claims  that  input  and  output  data  can  be  easily  changed 
to  arbitrary  units  (Ref.  1,  IV-C-4),  instructions  are  not  provided  how  to 
implement  such  changes,  nor  do  such  changes,  if  implemented,  affect  the 
Input  information  for  the  equations  of  state.  The  internal  calculations  are 
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-■ 


done  In  prescribed  units  which  were  also  adopted  for  the  U612  subroutine 
and  documented  In  the  comments  to  U6i2.  The  information  about  the  units  of 
all  quantities  is  Important  because  some  of  then 
(3<ln4)/3n,,  ar,/3n  ,  e,  e'  )  are  not  dimensionless,  contrary  to  statements 
in  the  TIGER  documentation.  Input  units  for  the  Lennard-Jones  (6-  2) 
potential  parameters  were  chosen  to  be  10*10  m  (Sngstrdms)  for  o  and  kelvins 
for  e/k.  A  change  of  these  units  would  require  minor  changes  In  the  U612 
subroutine  and  would  not  affect  the  rest  of  the  TIGER  program. 


In 

permits 

routines 


conclusion,  the  addition  of  the  two  subroutines  U612  and  UBS 
one  to  use  a  truncated  virial  equation  of  state  and,  because  the 
ace  carefully  commented,  facilitates  the  addition  of  further 


gaseous  state  equations. 
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APPENDIX  A 


LISTING  OF  TEE  SUBROUTINE  LJ612 


SUBROUTINE  LJ612 ILANE 1 


TO  COMPUTE  THE  OUTPUT  POP  STaTEG  USING  a  TRUNCATED 
VIRI AL  EQUATION  OF  STATE  IN  ahICh  THE  SECOnO  COEFFICIENT 
IS  9ASE0  ON  LENNARO-JONES  6-12  POTENTIAL  AND 
THE  THIRD  COEFFICIENT  IS  COMPUTED  USING  A 
SIMPLIFIED  RIGID  SPHERE  FORMULA 

INPUT  AND  OUTPUT  DEPEND  ON  The  VALUE  OF  LANE 

input  variables 

INPUT  FOR  LANES  1  THROUGH  » 

LO  A  NUMBER  OF  OUTPUT  UNIT  (PRINTER) 

NSG  A  number  of  GASEOUS  CONSTITUENTS 

RHO  *  mass  of  mixture  per  VOLUME  of  GAS  IN  G/CM*A3 

THIS  IS  CALLED  Rno-HAT  IN  TEXT*  SEE  I-B-I6. 
s*  3  sum  of  gaseous  mule  numbers  xid*  in  moles. 

T  A  GAS  TEMPERATURE  In  KELVINS 

JH  A  REFERENCE  mass  IN  grams  IM-xERO  In  I-B-16) 

X(U  A  mole  number  of  The  i-tm  gaseous  constituent 

MITH  1a1(...<  NSG*  IN  MOLES 

INPUT  FOR  LANES  S  AND  6 


INPUT  FOR  LANE  7 

IOEAL  *  NUMBER  OF  THIS  EQUATION  OF  STATE 
VALID  a  VALUE  OF  THE  SECOND  FIELD  OF  THE  "GEOS” 

INPUT  CARO 

INPUT  FOP  LANE  H 

VALID  A  VALUES  OF  THE  FIcLDS  2*3**  IFOR  Iai*2*3) 

OF  THE  "SET"  INPUT  CARO 

Input  for  lane  9 

VALID  a  VALUES  OF  THE  FIELOS  2  through  S  tl.l  THROUGH  *1 
OF  THE  "STG"  INPUT  CARO 


INPUT  FOR  LANE  10 

LO  a  NUMBER  OF  OUTPUT  UNIT  (PRINTER) 

PRNNTC  a  PRINT  I01CAT0R 

VALID  a  SIGMA  OF  A  CONSTITUENT.  IN  ANGSTROMS 
VALI2)  A  EPSILON/K  OF  A  CONSTITUENT*  IN  KELVINS 

INPUT  FOR  LANE  11 


a  INDEX  OF  »  CONSTITUENT 
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c 

c 

c 

c 


TYP(l)  a  SIGMA  OF  CONSTITUENT  NOALF  IN  ANGSTROMS 
7YPI2)  =  EPSILON/K  OF  CONSTITUENT  NGAlF  IN  KELVINS 

INPUT  FOP  LANE  12 

CO  *  NUMBER  OF  OUTPUT  UNIT  (PRINTER) 

NOALF  a  INOE*  NUMBER  OF  •  CONSTITUENT 

PRNNTC  =  PRINT  INDICATOR 

INPUT  FOR  LANE  13 

NOALF  =  INOEA  NUMBER  OF  A  CONSTITUtNT 
INPUT  FOR  LANE  1* 

LO  =  NUMBER  OF  OUTPUT  UNIT  (PRINTER) 


OUTPUT  VARIABLES 

output  for  lane  1 

CPHI  =  IMPERFECTION  FACTOR  PmI  In  THE  EQUATION  of  STATE 

PHIRHO  =  PARTIAL  DERIVATIVE  OF  LOG(CPHl)  mITM 
RESPECT  TO  LOGIfiMUl 

PHIT  a  PARTIAL  OERIVATlvt  OF  LOG(CPHl)  «ITH 

RESPECT  TO  LOG ! T ) 

OUTPUT  FOR  LANE  2 

INCLUDES  LANE  1  OUTPUTS  PLUS 

EPS  »  EnERuV  CONTENT  UF  GAS  CUE  TO  IMPERFECTION 

IN  MOLES.  (SEE  Il-C-11). 

EPSPT  =  PARTIAL  OERIVATIyE  OF  (T*EPSI  MlTM 

RESPECT  TO  T.  IN  MOLES. 

GAMRHO(I)  »  PARTIAL  DERIVATIVE  OF  GAMMAd)  .ITM 
RESPECT  TO  LOG(RmU) .  Ia1,....NSG. 

GAMT(I)  A  PARTIAL  DERIVATIVE  OF  GAMMAd)  »ITM 
RESPECT  TO  LOGIT).  . . . 

PMIN(I)  A  PARTIAL  DERIVATIVE  OF  LOG(CPMI)  WITH 

RESPECT  TO  (II)i  IN  I/MOLES.  1=1.. ...NSG. 

OUTPUT  FOR  LANE  3 

GAMMAd)  =  LOGARITHM  OF  THE  ACTIVITY  COEFFICIENT  OF 

THE  1-TM  CONSTITutNT.  1=1.. ...NSG.  (SEE  II-C-ll) 

OUTPUT  FOB  LANE  * 

Od.K)  A  PARTIAL  OEPIVATWt  OF  GAMMA(I)  .ITM  RESPECT 

TO  AIK)  IN  1/MOLES.  I.K=l*...*NSG.  (II-C-ll) 

GAMMAd)  A  LOGARITHM  OF  THE  ••CTIVITY  COEFFICIENT  OF 
THE  I-TM  CONSTITUENT.  1=1.. ..MSG. 

OUTPUT  FOR  LANES  S.B  ANO  11 

NONE 


LJ612 

59 

LJ612 

60 

LJ612 

61 

LJ612 

62 

LJ612 

63 

LJ612 

64 

LJ612 

65 

LJ612 

66 

LJ612 

67 

LJ612 

68 

Lj612 

69 

LJ612 

70 

LJ612 

71 

LJ612 

72 

LJ612 

73 

LJbl  2 

74 

LJ612 

75 

LJ612 

76 

LJ612 

77 

LJ612 

78 

LJ612 

79 

LJ612 

30 

LJ612 

61 

LJ612 

82 

LJ612 

83 

LJ612 

84 

LJ612 

85 

LJ612 

56 

LJ612 

87 

LJ612 

88 

LJ612 

89 

LJ612 

90 

LJ612 

91 

LJ6a2 

92 

LJ612 

93 

LJ612 

94 

LJ612 

95 

LJ612 

96 

LJ612 

97 

Lv612 

98 

LJ612 

99 

LJ612 

100 

LJ612 

101 

LJ612 

102 

LJ612 

103 

LJ6I2 

104 

LJ612 

105 

LJ612 

106 

LJ612 

107 

LJ612 

108 

LJ612 

109 

LJ612 

no 

LJ612 

111 

LJ612 

112 

LJ612 

113 

LJ612 

11* 

LJbl  2 

115 
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OUTPUT  FOR  LANE  6 

flOCTS  -  NUHBtR  OF  PARAMETERS  THAT  DEFINE  THE 

EQUATION  OF  StAT»-.  FOR  EACH  CONSTITJE*?  • 

FOR  THIS  EQUATION  NuCTS*2t  tM£  T-C 
PARAMfcTERS  ARE  LEnnAkO-J?** S  POTENTIAL 
PARAMETERS  SIGMA  mND  EPSILON/K.  SEE  LANE  II# 

OUTPUT  FOR  LANE  7 

KGEOS  =  NUMBER  Of  THIS  EQUATION  OF  STATE 
OUTPUT  FOR  LANE  9 

TYP(l)  -  AN  INDICATOR  FOR  SUBROUTINE  LIBRARY 

ZNAME  =  THE  value  OF  the  name  OF  A  CONSTITUENT 

VAL(3>  -  CEFAU.T  VALUE  OF  SIGMA 

VAL(A)  a  DEFAULT  VALUE  OF  tPSILOfl/K 

OUTPUT  FOR  LANES  10*12  ANO  1 * 

PRINTEO  OUTPUT  ON  UNIT  LO. 

OUTPU7  FOR  LANE  13 

TVPCl1'  =  SIGHA(NOALF)  IN  ANGSTROMS 
TYP (2)  =  EPSON  INQALF)*£PSIL0N/K  IN  KELVINS 

(KsRULTZMANN  CONSTANT) 

OUTPUT  FOR  LANE  GREATER  THAN  1* 

ERROR  PRINT  ON  UNIT  LO  ANO  STOP 


LOCAL  VARIABLES 


ALJ61 

ez 

CKA 

CZ 

EPSAB 

EPSOM  30) 

MISTAK 

NN 

RhOH 

SB 

SBKA(30I 

SBP 


3  CODE  VALUE  OF  THE  ALPHANUMERIC  NAME  OF 

THIS  ROUTINE  "LJ612".  THE  VALUE  IS  262*060102.0 
3  FACTOR  8-ZERO  OF  ThE  SECOND  VIRIAL  COEFFICIENT 
IN  ICm/ANGSTR0N)**3/MQLE.  SEE  LANE  S# 
a  THIRO  VIRIAL  COEFFICIENT  CF  CONSTITUENT  KA 
s  FACTOR  C-ZERO  OF  THE  THIRO  VIRIAL  COEFFICIENT 
IN  (CH/ ANGSTROM) **6/mOL£**2*  SEE  LANE  S# 

=  EPSILON/K  OF  THE  LENNARD-JONES  POTENTIAL 
sETwEEN  CONSTITUENTS  KA  AND  KB*  IN  KELVINS* 

3  awkAY  OF  EPSILON/«v  FOR  UP  TO  30  CONSTITUENTS 
s  ERROR  INDICATOR  IN  SUBROUTINE  LJttS  ARGUMENTS 
3  TEMPORARY  STORAGE  FOR  ERROR  PRINT 
3  RHO/«H  IN  l./CM*»j 
s  WEIGHTED  SUM  OF  o-STARS  IN 
MOl£S**2*AnGSTRQmS**3 
r  *>ARTIAL  WEIGHTED  SUMS  CF  9-STARS  IN 
MOLES* ANGSTROMS»*3 

s  WEIGHTED  SUM  OF  T-STAR  TIMES  B-STAR  DERIVATIVE 
IN  M0LES**c*ANGST°CM**3 


LJb'.Z 

116 

LJ612 

117 

LJol2 

118 

LJ612 

119 

LJ612 

120 

LJ612 

121 

LJol2 

122 

LJ612 

123 

LJ612 

12* 

LJ612 

125 

LJ6U 

126 

LJ612 

127 

LJ612 

12S 

LJ612 

129 

LJ612 

130 

LJ6I2 

131 

LJ612 

132 

LU612 

133 

LJ612 

13* 

LJ612 

135 

LJ612 

136 

LJ612 

137 

LJ612 

136 

LJ612 

139 

LJ612 

1*0 

LJ612 

1*1 

LJ612 

1*2 

LJ612 

1*3 

LJ612 

14* 

LJ612 

1*5 

LJ612 

1*6 

LJ612 

1*7 

LJ612 

1*8 

LJ6I2 

1*9 

LJ612 

ISO 

LJ6I2 

151 

LJ612 

152 

LJ612 

153 

LJ612 

IS* 

LJ6I2 

155 

LJ612 

156 

LJ612 

157 

LJ612 

158 

LJ612 

159 

LJ612 

160 

LJ612 

161 

LJ612 

162 

LJS12 

163 

LJ612 

16* 

LJ612 

165 

LJ612 

166 

LJSI2 

167 

LJ612 

168 

LJS12 

169 

LJ6I2 

170 

LJ612 

171 

LJ512 

172 

25 


c 

S8PKA<30!  =  PARTIAL  WEIGHTED  SUMS  IN  SBP 

LJ612 

173 

L 

HOLES* ANGSTROmS*"3 

LJ612 

174 

C 

S82P  =  wEIGHTEO  SUM  OP  <T-STaRI*"Z*(8-STaR-2-orIpEI 

LJ612 

175 

c 

IN  MULES** 2* ANGSTROM** 3 

LJ612 

176 

c 

SC  *  SEIGmTEO  SUM  OS  C-STARS.  mOlES*ANGSTROMS**6 

LJ612 

177 

c 

SIGA8  *  SIGMA  OP  THE  LENJ.ARO-JONES  6-12  POTENTIAL 

LJ612 

178 

c 

BETaEEN  CONSTITUENTS  KA  AnO  Kt>.  angstroms 

LJ612 

179 

c 

SIGA83  =  SIGAo**3.  ANGSTR0MS**3 

LJ612 

180 

c 

SIGPCT  =  REOUCTION  PaCTOR  UP  SIGMA  POR  ThE  RIGIO  SPHERE 

LJ612 

181 

c 

APPROXIMATION  IN  THE  TMIRO  VIRIAL  COEPPICIENT. 

LJ612 

182 

c 

OEFAULT  IS  0.81 

LJ612 

193 

c 

SIGMAC30)  *  -ELL  LOCATION  OP  lENNARO-JONES  6-12  POTENTIAL 

LJ612 

134 

c 

(IN  ANGSTROMS!  Puk  UP  TO  30  CONSTITUENTS 

LJ612 

185 

c 

TS  -  T-STAK.  ARGUMENT  IN  8-STAR  TABLES 

LJ612 

186 

c 

XA8  »  X<KAI*X(K9!  IN  M0lES**2 

LJ612 

187 

c 

LJ612 

188 

c 

LJbl2 

199 

c 

LJ612 

190 

c 

INPUT  NEEOEO  TO  ACTIVATE  THIS  EuuATION  OP  STATE 

LJ612 

191 

c 

LJ612 

192 

c 

LJ612 

193 

c 

1 

GEOS-INSTRUCTION  SPECIFYING  GASEOUS  EQUATION  OF  STATE 

LJ612 

194 

c 

CONSISTS  OP  THE  FOLLOWING  INPUT 

LJ612 

195 

c 

LJ612 

196 

c 

GE0S.LJ612 

LJ612 

197 

c 

LJ612 

198 

c 

2 

SET-INSTPUCTION  IF  THE  FACTOR  0.81  FOR  REDUCTION  OP 

LJ612 

199 

c 

SIGMA  IN  THE  RIGIO  SPHERE  APPROXIHATION  OP  THE  THlRO 

LJ612 

200 

c 

VIRIAL  COEFFICIENT  SHOULD  8E  CHANGEO.  SEE  LANES  S  AND  8. 

LJ612 

201 

c 

THE  INPUT  IS  AS  FOLLOWS 

LJ612 

202 

c 

LJ612 

203 

c 

SET.LJ612.SFACT.R 

LJ612 

20«* 

c 

LJ612 

205 

c 

WHERE  R  IS  THE  REPLACEMENT  VALUE  FOR  0.81 

LJ612 

206 

c 

LJ612 

207 

c 

3 

STO- INSTRUCT  ions  TO  SPECIFY  the  TWO  PARAMETERS  OF 

LJ612 

206 

c 

THE  LENNARO-JONES  6-12  POTENTIAL  FOR  EACH 

LJ612 

209 

c 

CONSTITUENT.  THE  INPUTS  ARE  AS  FOLLOWS 

LJ612 

210 

c 

LJ612 

211 

c 

STG.LJ612. NAME. SIGMA. EOK 

LJbl2 

212 

c 

LJ612 

213 

c 

WHERE  -NAME"  IS  The  ALPHANUMERIC  NAME  OF  THE  CONSTITUENT. 

LJ612 

214 

c 

"SIGMA"  IS  THE  SIGma  VALUE  In  ANGSTROMS  (1.0E-10  METRES! 

LJ612 

215 

c 

ANO  "EOK"  IS  the  value  OF  EPSILON/K  IN  KELVINS. 

LJ612 

216 

c 

K  8EING  The  Boltzmann  constant  1.380662E-23  J/K 

LJ612 

217 

c 

LJ612 

219 

c 

LJ612 

219 

c 

At VARS  CEIMINS  FECIT  29  APRIL  1983, 

LJ612 

22G 

c 

LJ612 

221 

c 

LJ612 

222 

COMMON  /  TIGER  / 

LJ612 

223 

X 

A(4l»10>*  AOEXP*  AJSilOl  •  AJ$J*  AL(40*10)* 

LJ612 

224 

X 

ALPHA*  AS  ( 1**«10)  *H  (40*12)  »  SETA*  *&(30)« 

LJ612 

225 

X 

oIG*  GJS(lU)*  RJSJt  CC(30*121 *CG(30*9)* 

LJ612 

226 

X 

CHII (30)  •  CHIJSU0I*CHIJSJ*  CHO (25) •  CL(1Q*9>* 

LJ612 

227 

X 

CN(iO)*  CNETAU0*12)*CPHl*  CPI  (30)  *  CPJS(10)> 

LJ612 

228 

X 

CPJSJ.  CS( 10*9) «  CV*  0(30*30) •  £«ED* 

LJ612 

229 

A 

ELH<10).  EOF.  EPETA1* 

EPETA2* 

EPETA3. 

LJ612 

230 

X 

EP£TA*».  EPS.  E°SPSI  * 

EPSPT. 

ETA ( 12) » 

LJ612 

231 

X 

ETETA1 (12) ,E0.  F (30) , 

FL. 

FLTA (50) 

LJ612 

232 

common  /  tigre  / 

LJ612 

233 

X 

FORM(SG,20) »r *£2(40) »G, 

GAMMA (30) 

«GAMPnO(30) . 

LJ612 

234 

X 

GAmT <30) ,  GETAU2).  h. 

HCON.HD* 

h£TA ( 12) . 

LJ612 

235 

X 

HETATl,  HETAT2*  MOF.HO. 

IUUG* 

1C  (82). 

LJ612 

23fc 

X 

IDEAL.  IEKROH.  IHIST* 

INOPT. 

INSTYP. 

LJ612 

237 

X 

ITC1L.  ITC2L*  ITC*L» 

I  TOP. 

IXCC(IO). 

LJ612 

238 

X 

IXMCJDS.  IXCTUO*.  JfMP* 

KGEOS* 

MILL* 

LJ612 

239 

X 

<IN0<1Q)»  KPPINT.  KX(30) • 

LI* 

LIB. 

Ljo12 

240 

X 

LINCT,  LO.LOPT,  MAXCHO. 

MAXFOP. 

MAXOPO. 

LJ612 

241 

X 

MAXPT*  MAXR£J.  MIXED* 

NCC* 

NCCP 

LJ612 

242 

COMMON  /  TYGEQ  / 

LJ612 

243 

X 

NCC1*  NCG.  NCT. 

NCT1. 

NCT2. 

LJ612 

244 

X 

NNIS.  NOalF.  NOCASE* 

NOCK). 

NOCTS. 

LJ612 

245 

A 

NOFLT*  nO^Omh.  NOGCCIO) 

•  NOGEOS* 

NGORO. 

LJ6I2 

246 

X 

NOREJ*  NOTOT*  NPAGE* 

NPTSV • 

NSC. 

LJ612 

247 

X 

NSG.  NSM*  NST* 

NTOC* 

ORU(2S>. 

LJ612 

248 

A 

P.  PCON.  PETAU2) 

♦  PETAT1* 

PETAT2. 

LJM2 

249 

A 

PmIJSCIOI.PHIJSJ*  PHIN(30> 

t  PHIPmO. 

PhIT. 

LJ612 

250 

X 

PSI  (30) «  PSU.  PSI22. 

PO. 

0(10). 

LJ612 

251 

X 

OA(IO) •  R.  R£J (SO) • 

RHO. 

RhOO 

LJ612 

252 

COMMON  /  TYGae  / 

LJ6I2 

253 

X 

Px*S*SO*  SENT  , scnTE, SG (30) * 

SQQ* 

STOCN(IO.U)* 

LJ612 

2S4 

X 

STOFl  'll)  .STOGnUO.II)  .STGMM11) « STONX  l 30 ) . ST OPC  110*11)* 

LJ612 

25S 

X 

SVOL*  SVOLE* SvP (200*2) ,SX,$0, 

T. 

TACON*  T£ma (9.20) * 

LJ612 

256 

X 

TInITL.  TMCON.  TSAVE* 

TYP(«2) • 

TO. 

LJ612 

2S7 

X 

umi  (3o> *  umjsuo)*  umjsj* 

V. 

VAL (82) * 

LJ612 

258 

X 

VCON*  v£TA (12) «  VETAT1* 

VETAT2* 

VO. 

LJ612 

259 

X 

M(10),  NCrll.  WCP* 

MM. 

MMU. 

LJ612 

260 

X 

*P.  MT(4Q),  X  (30)  • 

JETJTl 170) .XETAT2 130) .XN. 

LJ612 

261 

X 

xnahE (40)  •XUM1*AUM2*Y(30) .ZC (82) *ZCC<80) .ZCx<80> ,  ZNAmC 

LJ612 

262 

c 

LJ6I2 

263 

COMMON  /  P»INTC  /  NGPRCN.PRNNTC 

LJ612 

264 

LOGICAL  PRNNTC 

LJ612 

265 

c 

LJ612 

266 

DIMENSION  EPSOKl3QI*$IGMA<3Q).Sd»vA(30),SBPiU(30> 

LJ612 

267 

c 

LJ612 

268 

DATA  AlJ612/262<*C60102.0/ 

LJ612 

269 

c 

*LJ6l2  IS  THE  PEPPEStNTATION  OF  T«t 

ALPHANUMERIC 

LJ612 

270 

c 

"LJ612**  IN  THIS  COOE.  (SEE  SU8«QUTJnE  CARQfiO*) 

LJ612 

271 

c 

LJ612 

272 

lFd.jNE.GT.  1*16010  1S00 

LJ612 

273 

GOTO!  100, 100*30  0*<*U  0*50  0*60  0*700 

.800*400*1000  «U00*  1200* 

LJ612 

274 

* 

1300*1400) *LAn£ 

LJ612 

275 

c 

LJ612 

276 

c 

LANES  I  ANO  2 

LJ612 

277 

c 

LJ612 

278 

100 

JO  104  KAsl«NSG 

LJ612 

279 

IF(SIGMAlKA)  •GT.O..ANO.EPSOMKA) 

«GT«0.)  GOTO  104 

LJ612 

280 

c 

LJ612 

2*1 

101 

SN=101 

LJ612 

282 

xRITE  (LO« 102)  NN.AA.SIGMAIKA)  ,KA 

*E°SOK(KA)* 

LANE 

LJ612 

283 

102  FORMAT  1 lHQ, 10X,31H£^MCP  STOP  IN 

EQUATION  OF 

STATE, 

LJ612 

284 

A 

2SM  POUTIkE  LJ612  AT  STATEMENT 

» I4.9h  BECAUSE./ • 

LJ612 

285 

8 

1H  » 10X«6mSIGmA( *I2«3H)  =.1PE12 

•S.SA.IOHANO  EPSOK1.I2. 

LJ612 

286 

27 


C  3H<  *.IPE12.S*5X.21MARE  NOT  BOTH  POSITIVE*/* 

0  IH  .10X.6HUNE  **131 
STOP 
C 

104  CONTINUE 
C 

S0  *  0, 

S8P*  0, 

SC  *  0. 

SB2P*  0. 

C 

00  145  KA*l*NSG 
C 

IF (LANE.NE.2)GOTO  105 
S8KA(KA>*0. 

S8PKA(KA)aO* 

C 

105  SC*SC»X(KA)«SIGKA(KA)*«6 

c 

DO  135  K8*1*NSG 
C 

SIGA8*(SIGmA(KA)«SIGMA(KB)}/2. 

SIGA83*SIGA5**3 

EPSA8*S0R  T <  C  < SIGMA ( K A) /SIG AB>  » < SI GMA  t K0 ) /SI GAB > ) **3 

*  *epsok<ka>«epsok man 

XA8*X(KA)*X<K8) 

TS=T/£PSA8 

c 

CALL  LJBSlLANE*TS*dS*T0SP*TT9SPPtMlSTAK) 

C  INTERPOLATE  IN  THE  8STAR  TABLES 
IF(HJSTAK,EQ.O)GOTO  125 
111  NN*111 

WRITE  (L0*U5lNN.MlSTAKtLANE*TS 

-RITE (LO. 116} KA.K8.T 

STOP 

115  FORHATC1H0.10X.31HERR0R  STOP  IN  EOUATION  CF  STATE. 

*  42H  ROUTINE  LJ612  AFTER  CALLING  INTERPOLATION. 

♦  27M  ROUTINE  LJ0S  AT  STATEMENT  *14./. 

*lh  ,10X.8hmISTAK  *.I2fSX«6HLAN£  =, I4.SX.4MTS  *.1PE12.5> 

116  FORMAT  (1M  .10X.4MICA  *.I3*5X.4HKe  *.I3«5X* 

♦  3MT  *.1PE12.5} 

C 

125  S8*S8 ♦ S I GA03#BS*X AB 

S8P=S8P*SIGAB3*T8SP*XA8 
IF (LANE. NE.2) GOTO  135 

C 

S3KA (KA) *S8KA (KA) *SIGAB3»8S*X (Ko) 
SBPKAIKAI*SBPKA(KA)»SIGAB3*TBSP*A(Kd) 
SB2P*SB2P*SIGA03*TTBSPP*XA8 
C 

135  CONTINUE 
C 

l-*5  CONTINUE 
C 

RhOMxRhO/WM 

S8*82*S8*RhOM/SX 

SC*CZ*SC*RM0M**2*SA 


LJ612 

287 

LJ612 

288 

LJ612 

289 

LJ612 

290 

LJ612 

291 

LJ612 

292 

LJ612 

293 

LJ612 

294 

LJ612 

295 

LJ612 

296 

LJ612 

297 

LJ612 

298 

LJ612 

299 

LJ612 

300 

LJ612 

301 

LJ612 

302 

LJ612 

303 

LJ612 

304 

LJ612 

305 

LJ612 

306 

LJ612 

307 

LJ612 

308 

LJ612 

309 

LJ612 

310 

LJ612 

311 

LJ612 

312 

LJ612 

313 

LJ612 

314 

LJ612 

315 

LJ612 

316 

LJ612 

317 

LJ612 

318 

LJ612 

319 

LJ6I2 

320 

LJ612 

321 

L  J612 

322 

LJ612 

323 

LJ612 

324 

LJ612 

325 

LJ612 

326 

LJ612 

327 

LJ612 

328 

LJ612 

329 

LJ612 

330 

LJ612 

331 

LJ612 

332 

LJ612 

333 

LJ612 

334 

LJ612 

335 

LJ612 

336 

LJ612 

337 

LJ612 

338 

LJ612 

339 

LJ612 

340 

LJ612 

341 

LJ612 

342 

LJ612 

343 

28 


m 


CPMI=S8*SC*1.0 

PhlRH0=(S8*2.*SC)/CP«I 

SBP=9Z*S3P 

PMlTs»HOM*S0P/(SX*CPMl) 

IP (LANE.EQ.U  GOTO  1900 
C 

C  WANE  2 
C 

200  £PS=-RHCM»SdP 

£3SPT=-Rhoh* (S9P#2. ♦9Z*S32°) 

c 

00  255  KA=1*NSG 
C 

CKAsC2«SXGHA«KA)*»t> 

GAMRHO(KA>=2.«SC*2-#KHO»*HZ*SbKA(KA)* 

*  (RhO**SX) ••2*CXA 

GAMT (KA) 52.»RHOM#BZ*SBPKA (KA) 

PHlN(KA»a{J-S3*SC)/SX*2.*RM0H*bZ*SBKAlKAI/SX 

•  .RH0H««2*SX*CKA)/CPtiI 
C 

255  CONTINUE 
GOTO  1900 
C 

C  LANE  3 
C 

300  OO  305  KA=1#NSG 

IF(SIGMA(KA).GT.l)..ANO.EPSOK(KM.GT.O.)  GOTO  305 

c 

30*  '/N=30* 

-RITE (LO* 102)  NN.KA.SIGMA (KA) .KA.EPSOK (KA) *LAN£ 

STOP 

C 

305  CONTINUE 
C 

SC-0 

OO  3*5  KA=1*NSG 

c 

SoKA(KA)=0. 

GAMMA (KA) =SIGKA (KA) **6 
SC*SC*GAMMA (KA) *X (KA) 

C 

00  335  KB=1*NSG 
C 

SIGA9= (SIGMA (KA) •5Ito«A(K8) )/2. 

SIGAB3*SIGAB»*3  „  ^ 

£PSA8=S0RT( ( ( SIGmA ( KA ) /SI 5AB) •( SIGMA (KB) /biGAB) )*«3 
♦’•EPSOK  (KA)  •EPSOK  (*<IJ) ) 

TS*T/EPSAB 

C 

CALL  L J0S (LANE* T**dS* T0SP*TT8SPP« HIST AK) 

C  INTEOPCLATE  IN  THE  3STAH  TABLES 
IF(mISTAK.£0.0)ooT0  325 

C 

311  NN=3H 

k*UTE  (LO*  US)  M*nISTak*lANE«TS 

•RITE (L0*116'KA*Kd*T 

STOP 


29 


C  STOP  IF  TS  IS  OUTSIDE  ThE  RANGE  OF  Tm£  8STAR  TABLES 

C 

325  S9KA (KA J  =SBKA (KA) »SIGA83*0S*X <KO) 

C 

335  CONTINUE 
C 

3<*S  CONTINUE 

r 

RMQM=RhO/-M 

00  355  KA=1.\SG 

C  GAKMA 1KAI  =  {  (GAMMA  (KA)  •(!  .5#SX*SC)  •SX»RMOM*CZ* 

♦  2.*SBKA«AI*BZ)*RrtOM 

C  GAMMA (KA)  ALREAOY  CONTAINS  C(KA) 

C 

355  CONTINUE 
C 

GOTO  1900 
C 

C  LANE  •* 

"  AOO  SC=0 

00  *15  KA=1*NSG 

C  IF(SIGMA(KA).GT.O*«ANO.EPSOK(KA).GT.O.)  GOTO  MO 

C 

fc09  sN^409 

KRITE  (LO*  102)  NM.KAtSIGMA  (KA)  »KA*£PS0K  (KA)  *L*Nt. 

STOP 

C  *10  SC=SC*X(KA) •SIGMA (KA)««6 

C 

4,15  CONTINUE 
C 

RHOMsRHO/hM 

C 

00  445  KArl*NSG 
C 

S9KA (KA) =0 
CKA=SIGKA (KA) +*b 
GAMMA (KA)sCKA 
C 

00  435  KB=l*NSG 

SIGABs(SIGMA(KA) ♦SIGMA(K8) )/2* 

EPSA8=S0RT  I  (  (SIGMA (KA)/SIGA8)#(MGmA(K8)/S1&AB)  )  * 
•  •EPSOK(KA)«EPSOK(K0)  I 
TS=T/EPSAP 

C  CALL  LJBS(LANE»TS*dS*TBSP*TT9SPM*M*STAK) 

C  INTERPOLATE  IN  THE  dSTAM  TABLES 
IF(MISTAK.EO*0)SOTO  «*25 
C 

4,21  nn=*21 

-KITE (LO*115)NN«m1STaK,LA,C»TS 
•RITE (LO* llG)KA*Ktf • T 


C  STOP  IF  TS  IS  OUTSIDE  T ♦€  INTERPOLATION  TABLE 

C  <*25  0CKA,KaJ=((SC*SA*(CKA*SIGNA(KB)**6n*CZ*RH0N 
♦  ♦2.*BZ*BS«SIGAB3>*ftHOH 
SBKA(KA) *SBKA(KA) ♦SIGAB3*BS*X (K6J 
C 

<*3S  CONTINUE 
C 

4*<*5  CONTINUE 


c 

DO  455  KA*1«NSG 

GAMMA  CKAJ a ( (GAMMA  CKA) #Q.S*SX*SC) •SX*CZ*RHOM» 
A  2,*SBKA(KAJ*BZ)*RrtOM 
455  CONTINUE 
C 

GOTO  1900 


C 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 


LANE  5 

500  BZ*U 261354 
SIGFCT*0.81 

CZ*BZ*#2*  (S./8. )  •SIt»FCT#*6 
GOTO  1900 

urne  s  is  called  from  tigera-io  to  establish  default  values  of 
constants  In  the  eouation  of  state 

the  SECOND  VIRIAL  COEFFICIENT  IS  B*3STAR*8Z*SIGNA**3. 

THE  OIMENSION  OF  82  IS  1CM/ANGSTR0NI ••3/NOLE. 

REQUIRING  SIGMA  IN  ANGSTROM*  RHO  IN  G/CN**3  AND 
HASS  «N  IN  GRAMS.  THE  VALUE  OF  C2  IS  ACCORDING  TO 
HIRSChfELOER*  CURTIS  AND  BIRO*  SECTION  III-5-A.  ANO  The 
THIRD  VIRIAL  COEFFICIENT  IS  C»CSTaR«C2*SIGHA**6. 

C2  MAY  BE  CHANGED  IN  LANE  B. 

LANE  6 

600  N0CTS=2 
GOTO  1900 

LANE  S  IS  CALLEO  from  T1GERA-10.  NOCTS*2  INDICATES 
THAT  EACH  CONSTITUENT  IS  CHARACTERI2ED  BY  TvO 
PARAMETERS  NHICH  SPECIFY  ITS  LENNAKD-JONES  6-12  POTENTIAL 


LANE  Y 

700  IF<VAL<U.EO.ALJ*121*GEOS»IOE*L 
GOTO  1900 

CALnEO  FROM  TIGERA-N20  THIS  INDICATES  IF  THE  6E0S-CARD 
CONTAINS  THE  NAME  OF  THIS  R0UTINEI-LJ612- 

LANE  B 

BOO  IF(VALII).NE.ALJ612100T0  1900 

IFIVALI21 .NE.33201S173A.01GOTO  1900 
SIGFCT=ABS1VAH3)I 
C2*9Z»»2,<S./8.1*S1gFCT**6 
GOTO  1900 

THIS  IS  CALLEO  FROM  TIGEPA-AAO  ANO  IT  CHANGES  THE 


L4612 

458 

LJ612 

459 

LJ612 

460 

LJ612 
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LJ612 

496 
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504 
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C  OEFAULT  FACTOR  0.81**6*0.2824  IN  THE  CONSTANT  C2. 

C  THE  CODE  OF  VAL<2)  IS  MSFACT"  WHICH  MUST  BE  ON 

C  THE  input  CARD  AFTER  MS£T.LJ612*H.  NOTICE  THAT  THE 

C  INPUT  IS  THE  REOUCTION  FACTOR  OF  SIGMA  ANO  IS  RAISEO  TO  6-TH  POWER 


900  IF (VALID .NE.AL J612)G0T0  1900 
TYP(1>=2.0 
ZNAM£sVAL(2) 

IF (VAL (3) .LE.O. )  VALO)  =  3.5 
IF(VAL(4) ,L£.0.>  WAL14)  a  300. 

GOTO  1900 

THIS  IS  A  CALL  FROM  SUBROUTINE  LIBRARY-400  (III-C-245) 
PROCESSING  A  STG-CARD.  AT  THIS  TIhE  VAL (2)  CONTAINS 
THE  NAME  OF  ThE  COMPONENT.  VAL(3)*SIGMA.  VAL (4)*£PSOK. 
SIGMA  AND  EPSOK  ARE  REPLACED  8T  OEFAULT  VALUES  IF  THEY 
ARE  NOT  POSITIVE. 


1000  IF CPRNNTC)  WRITE (LQ. 1 015)  VAL(1)*VAL(2) 

1015  FORMAT (1h  ,7X.5HLJ612.3X.6HSIGMA3#IPE12.5.10H  ANGSTROMS* 
♦  15H,  EPSIL0N/K*.1PE12.S.8H  KELVINS) 

GOTO  1900 

C  THIS  LANE  IS  CALLED  FROM  LIBRARY-940  (SEE  III-C-2S0) 

C  WHICH  HAS  JUST  STORED  SIGMA  ANO  EPSILON/K  IN  VALID  ANO 
C  VAL (2)  FOR  PRINTING 


1100  SIGMA(NOALF>*TYP(ll 

IF(SIGMA (NOALF) .LE.O. )  SIGMA(NOALF) *3.5 
EPS0K(N0ALFI*TYP(2) 

IF (EPSOK (NOALF) .LE.O.)  EPSOK (NCAlF) *300. 

GOTO  1900 

C  LANE  11  IS  CALLED  FROM  SUBROUTINE  COmPOS-530  (III-C-77) 
C  AND  5U9R0UTINE  R0RDE«-90  (III-C-312).  WHICH  HAVE 
C  STORED  THE  PROPER  VALUES  IN  NOALF  AND  TYP.  BY  THIS 
C  CALL  THE  ARRAYS  SIGMA  ANO  EPSOK  AR£  FILLED  UP  AND 
C  REORDERED*  RESPECTIVELY. 

C  SIGMA  ANO  EPSOK  ARE  REPLACEO  BY  THEIR  OEFAULT  VALUES  IF 
C  THE  VALUES  IN  TYP  AR£  NOT  POSITIVE. 


1200  IF(PRnnTC)  W«ITE(LO*121S)  SIGMA (NOALF) .EPSOK (NOALF) 

1215  FORMAT (IH  *13X*5HLJ612.3X«6hSlGMA*.lPE12.5*10H  ANGSTROMS* 

♦  15H,  EPSIL0N/K*.1PE12.5*8H  KELVINS) 

IF ( PRNNTC. ANO.S IGhm (NOALF ) .EQ.3.5. ANO .EPSOK (NOALF } .E0.300 . ) 
A  WRITE (LO. 1216) 

1216  FORMAT (1H*.90X*16H(DEFAULT  VALUES)) 

GOTO  1900 

C  LANE  12  IS  CALLED  FROM  COmPOS-860  (III-C-BO)  WHICH 
C  PRINTS  A  LIST  OF  GASEOUS  CONSTITUENTS 


c 

LJ612 

S72 

1300  TYPC1)*SIGmA(NOALF> 

LJ612 

S73 

TYP(2)*EPSOK(NOALF> 

LJ612 

S74 

GOTO  1900 

LJ612 

575 

C  THIS  LANE  IS  CAUEO  FROM  ROROER-40  UII-C-311)  MHICH 

LJ612 

576 

C  REORDERS  THE  LIST  OF  CONSTITUENTS*  LANE  13  IS 

LJ612 

577 

C  THE  INVERSE  OF  LANE  11* 

LJ612 

578 

C 

LJ612 

579 

C  LANE  14 

LJ612 

580 

C 

LJ612 

531 

1400  «RITE(L0#141S)SIGFCT 

LJ612 

582 

1415  FORMAT ( 1H0 • 30X • 35HTRUNCATED  VIRIAL  EQUATION  OF  STATE  • 

LJ612 

583 

♦  18H( SUBROUTINE  LJ612J./.1M  *10X* 

LJ612 

584 

•  S3HSECOND  VIRIAL  COEFFICIENT  COMPUTED  FROM  LENNAR0- JONES* 

LJ612 

585 

*  1SH  6-12  potential*/*im  iox* 

LJ612 

586 

♦  52HTHIRD  VIRIAL  COEFFICIENT  COMPUTED  USING  A  SIMPLIFIED. 

LJ612 

587 

♦  37H  RIGIO  SPHERE  FORMULA  MITH  R  *  SIGMA*.1P£12.5) 

LJ612 

S88 

GOTO  1900 

LJ612 

589 

C  LANE  14  IS  CALLED  FROM  PRINNT-SO  CI1I-C-297) REQUIRING 

LJ612 

590 

C  TO  PRINT  GENERAL  INFORMATION  ABOUT  THE  EQUATION 

LJ612 

S91 

C  OF  STATE 

LJ612 

592 

C 

LJ612 

593 

C  LANE  GREATER  THAN  14 

LJ612 

594 

C 

LJ612 

595 

1500  1RITECL0.15ISILANE 

LJ612 

596 

STOP 

LJ612 

597 

ISIS  FORMAT! 1H0.10X.32HSTOP  0Y  SUBROUTINE  LJ612  BECAUSE* 

LJ612 

598 

♦  2SH  IT  WAS  CALLED  MITH  LAN£«*l5»/*lH  *10X* 

LJ612 

599 

•  26HLANE  SHOULD  BE  LESS  THAN  15.) 

LJ612 

600 

C  ERROR  STOP  MITH  INVALIO  VALUE  OF  LANE 

LJ612 

601 

C 

LJ612 

602 

1900  RETURN 

LJ612 

603 

END 

LJ612 

604 
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LISTING  OF  THE  SUBROUTINE  UBS 


APPENDIX  B 


LISTING  OF  THE  SUBROUTINE  UBS 
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c 

c 

c 
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c 

c 
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SUBROUTINE  I _ IBS lLANt»TS«8St T9!»P*7T8SPPf HlSTAKl 

THIS  ROUTINE  PROVIDES  INTERPOLATED  VALUES  OP  BSTAR 
MO  ITS  DERIVATIVES.  IT  IS  CALLED  PY  THE 
T«U*<CAT£3  VjaiAL  SUdROUTlNE  LJ612.  dASEO  ON 
LENNARO-JCJ.ES  (6-12)  POTENTIAL. 

LANE  IWICATES  -MAT  SHOULD  BE  INTERPOLATED 

=  1  COMPUTE  dSTAR  and  FIRST  DERIVATIVE 
=  2  COMPUTE  dSTAR  ANO  FIRST  TWO  DERIVATIVES 
a  3  COMPUTE  BSTAR 
=  A  COMPUTE  6STAR 
TS  =  ARGUMENT  TSTAR 

The  FOLLOWING  WILL  t)£  COMPUTEO  Bv  this  ROUTINE 

as  *  esTAR 

T9SP  =  TSTAR*IFIRST  DERIVATIVE  OF  BSTAR) 

TTBSPP  =  T ST Afl»*2* (SECOND  OERIVAllVE  OF  6STAR) 

MISTAK  =  ERROR  INDICATOR, 

s  0  If  NO  ERROR 

*  -1  IF  TSTAR  IS  LESS  Than  0.01 
a  2  IF  LANE  IS  OUTSIOE  INTERVAL  U»A) 

the  interpolation  is  oone  in  hirscmfeloer-cuRtis-biro 

tables  i-b  using  twu-point  herhite  formulas  if  o.?s.le.ts.le.S.o  . 

outside  this  range  series  expansions  are  useo. 

AI VARS  CELMINS  FECIT  13  APRIL  1983. 

DIMENSION  TRLJTa<82)*BSTLJ(82)*aSTLJl(82)«BSTLJ2(82) 

ThE  FOLLOWING  ARE  HIRSCHFELOER-CURTIS-BIRD  TABLES  I-B 

DATA  (TRLJTBH)  .1*1*82)  /.30*  .35.  .40.  .45*  ,50.  ,5S.  .60*  .65.  . 
170.  .75.  .SO.  .AS.  .90.  .95.  l.uO.  1.0S.  1.10*  1.15*  1.20*  1.25.  1 

2.30.  1,35.  1.40.  1.45.  1.50.  l.SS.  1.60.  1.65.  1.70.  1.75.  1.80.  1 

3.85.  1,90.  1,95.  2.00.  2*10.  2,20.  2,20.  2.40.  2.50*  2.6*  2.7.  2.8 

4.  2.9.  3.0.  3.1*  3.2.  3.3.  3.4.  3,5.  3,6.  3.7.  3.8.  3.9*  4,0.  4,1. 

S  4,2,  4,3,  4, A.  4«S.  4,6.  4,7.  4,8.  <*,9.  5.0.  6,0.  7,0.  8,0.  9,0* 
810,0.  20,0.  30,0*  40,0.  50,0.  60.0,  70,0.  80.0.  90,0.  100,0.  200,0 
7.  300,0.  400,0/ 

DATA  (BSTLJ(I) .1=1.82)  /-27.8ttOSolO«  -18,7548950*  -13.79883S0.  -10 
1.7549750*  -8,7202050*  -7,2740858*  -6.1979708*  -5.36S1918.  -4.710«3 
270.  -4,1759283.  -3,7342254.  -3,3631193.  -3.0471143,  -2.7749102,  -2 
3.S380B14,  -2.3302208.  -2. 14^3742,  *1,9825492.  *1.8359492,  -1.70377 
*»84,  *1.5841047,  -1,4752571.  -1.3758479,  -1.2847160.  -1.2008832.  -I 
5.1235183.  -1.0519115*  -.98S4S34,  -.9236164,  -.8659428.  -.8120333. 

6- . 7615373,  -.7141473.  -.6695903,  -.6276254.  -.5506331.  -.4817100. 

7- , 4196776,  -.3635757,  -.3126134,  -.2661334,  -,2235863.  -.184S073. 

8- . 1485022.  -.1152339.  -.0844124,  -.0557870,  -.0291400,  -.0042809, 
9.0189568.  ,0407201.  ,0611388.  *0603279,  , 098390 t.  .1154169,  .131*9 
A02.  .1466837.  .1610638.  .1746904,  .1876177,  , 1998951.  .2115673.  .2 
9226751,  .2332SS8*  .2433435.  .3229044,  ,3760885.  .4134340.  .4405978 
C,  .4608753,  .52S3742.  .S2692SS.  .S16S?S0.  .5083614.  ,4962126.  .488 
06507,  ,4797901.  .4716150.  .4640695,  .4114317,  .3801279.  .3583S12/ 

DATA  18STLJ1 tl),l*l,62)  /76.6072S60.  45.2477130.  30.2670800.  21.98 
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194020,  16.9236900.  13.S821S60*  11.2406490*  9.54S5Q96,  0.257114S*  7 
2.2540135.  6.4541400.  5.8034061.  5.2649104*  4.8127607*  4.4202616*  4 

3.09766S9.  3.8106421*  3.559292S*  3.3374893*  3.1404074.  2.9642040*  2 

4.6057826.  2.6626207*  2.5326459.  2.4141403*  2.3056683*  2.206021S.  2 

5.1141772*  2.0292621.  1.9505276.  1.0773207.  1.8091057.  1.7453722.  1 

6.6b570lo.  1.6297207.  1.5275444*  1.4366294.  1.3552180.  1.2819016.  1 

7.2155320.  1.1S51691.  1.1000353*  1.0494802*  1.0029572*  .9600031*  .9 

8202229.  .8832774.  .84 88746.  .8167606*  .706714S*  .7585430*  .7320758 
9,  .7071630.  .6836715.  .6614830*  .6404922*  .6206045.  .6017352.  .563 
40082.  .5667545.  .S50SU8.  .53S0237.  .5202387.  .S061101.  .4Q25951. 
8.3839722*  .3082566.  .2524801*  .8097011.  .1758670.  .0286638.  -.0174 
C929.  -.0393115.  -.0516478.  -.0S93621.  -.0645039.  -.0680819.  -.0706 
0470.  -.0725244,  -.0775400,  -.0765245.  -.0747534/ 

OATA  (BSTLj2(I) .1=1*82)  /-356.87o790.  -189.465360,  -116.366040.  -7 
18.8779500.  -57.3395200.  -43.8024500.  -34.9186900,  -28.6405000,  -24 
2.0626600,  -20.41311U0.  -17.94l9uOO*  -IS. 8254600.  -14.US5700,  -12. 
37108100,  -11.5398500.  -10.5513300,  -9.7074400,  -8.9798500,  -0.3470 
4000.  -7.7921700.  -7.3022700.  -6.8669200,  -6.^777700,  -6.1280500.  - 
55.dl22500,  -5.5257800.  -5.2648S00,  -5. 0262000.  -4.8073800,  -4.6058 
6700,  -4.4198000,  -4.2475000,  -4.0075300*  -3.9386300.  -3-7997200,  - 
73.5481400.  -3.3264700.  -3.1297400,  -2.9540100,  -2.7961400,  -2.6535 
6500.  -2.5241600,  -2.4062300.  -2.2983100.  -2.1992000,  -2.1070500,  - 
92.0234000.  -1.9451100,  -1.872310(1.  -1.8044700,  -1.7410000.  -1.6017 
A400.  -1.6260500.  -1.5737100*  -1.5244100*  -1.4778900*  -1.4339400*  - 
91.3923400*  -1.3529100,  -1.3154800,  -1.2799100,  -1.2460600,  -1.2138 
C100.  -1.1830500.  -1.1536700,  -.9193930,  -.7579300,  -.6396790.  -.54 
097920*  -.470779U*  -.1704030*  -.0720120*  -.0241090.  .0039270*  .0221 
£470.  .0340170.  .0440560*  .Q510310*  . 0564*10.  .0772960,  .0813970.  . 
F0820550/ 


IFCLANE.6E.1.AnD.LAN£.lE.4)GOTO  15 

HISTAK=2 

RETURN 

15  MlSTAKsO 

IFtTS.GE.TRLJTBUOl  .AnO.TS.LT. TM.JT8165J)  GOTO  2S 

IFUS.GE.0.01)  GOTO  36 

MlSTAKa-l 

RETURN 

NEXT  FIND  SLOT  FOR  INTERPOLATION 
25  00  35  KA=11,65 

IF  ITS.LE.TflLJTfl  <I<AIJ<50T0  45 

35  CONTINUE 

36  CONTINUE 

START  COMPUTING  FUNCTIONS  9Y  SERIES  EXPANSION 
Al*1.22S  417  702  7  •  SORT <2. ) •TS**C-0.25) 
A2=-3.t>25  609  900  •  SORT<2.1»0.s»TS«#<-0.75) 

A1  AND  A2  ARE  FIRST  T*0  TERMS  OF  SERIES  FOR  ttSTAR 
AKs-1 
8S=AI*A2 

T8SPs<-Al-A2*3.J/4. 

TT&SPP* I A1*S«« A2*2l . ) /16. 


3?  AK*AK.l* 


COMPUTE  NEXT  TNO  TERNS  OF  SERIES 
KVIT*1 

STOP  SERIES  COMPUTATION  WHEN  KVIT.NE.O 

AlsAl*  (U.O^AK-l. 01/ 12.0*  AX‘1.0)  1/(12. 0«AK»2.  01  »TS> 

A2*A2*I I*.0»AK«1.0)/(2.0»AK*2.0I)/I (2.0*AK*3.0) »TS> 

08>A1«A2 

es*as*o0 

IFU8SI08).OT.I.0£-12)  KVIT»0 
IF(UN£.GT.2I  GOTO  39 
COMPUTE  ONLY  8STAR  IF  LANE  *  3  0»  * 

08P«-(AK*1.2S)«A1*<AK.1.TS)»A2 

T8SP»TBSP*DBP 

1FIA8S108P) .GT.1.0E-12I  KVIT*0 
IFILANE.EO.lt  GOTO  39 

COMPUTE  ALSO  SECOND  DERIVATIVE  IE  LANE  »  2 

0BPP=><AK*l.2S)«<AK*2.2S)»Al»(AK*l.?S)»(AK»2.7S> *»2 

TTBSPP«TT8SPP»DBPP 

1F(A8SI08PP).GT.1.0E-12)  KVIT*0 

39  IF(KVlT.NE.O)  GOTO  *0 
BRANCH  IF  COMPUTATION  CONVERGED 

IFIAK.LE.30S.)  GOTO  37 

COMPUTE  AT  MOST  30S  TERMS  OF  SERIES  (NEEDED  FOR  TS*0.0l> 

40  RETURN 

ENTER  4S  FROM  35  AND  INTERPOLATE 
45  Tl»TRLJT8(KA-n 
T2«TRLJTB(KA1 
0ELTST*T2-TI 
XI* (TS-T1 1 /OELTST 
X2»(TS-T2)/0ELTST 
FI=(1.-2.*X1)«X2*«2 
F2m*Im*2m(I.-2.»X2) 

01=>X1»X2«42 

02aXiM2*x2 

FltFZt 01.02  ARE  HERHlTE  INTERPOLATIOH  FACTORS 
BS-F1*BSTLJ (KA-1>  *F2*0STLJ (KA) * 

♦  (Dl*BSTLJl(KA-l>/Tl*D2*BSTLJH*A)/T2)*0£LTST 
IF (LANE. 6T.2) RETURN 
RETURN  IF  LAN£33  0*  LANE«» 

ELSE  COMPUTE  FIRST  DERIVATIVE 

TSSP=F1*°STLJ1  '."r.A-I!  ♦F2*eSTLJM"'A>  *  DELIST* 

A  (01* I8STLJI IKA-l ) ♦BSTLJ2 (KA-l> I /T1 *02* (BSTlJI tKA) ♦BSTLJ2 (KA) ) 
IF ILANE.EQ. II RETURN 

IF  LANE*2?  ALSO  COMPUTE  SECOND  DERIVATIVE 

TT8SPP*<8STLJl I KA-ll -BSTLJI (KAI  >*6«*Xl**2*TS/0ELTST  * 

A  HBSTLJlCKA.n-faSTLj2IKA.l)l/Tn*TS.*2*0.*XI-U)  ♦ 

B  ((BSTLJI (KA) ♦BSTLJ2(KA) ) /T2I *TS*Xl*(3*#x2*l«)  -  TBSP 

RETURN 
EM) 
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SAMPLE  OUTPUT 


1 

! 


APPENDIX  C 
SAMPLE  OUTPUT 


lira: 

UJUul 

030 


a  x  x 
oo  o 


TXT 

O 


arTTattaaaxatxaerxzftzaaa: 

UitJUUUUlilUUlWWUUIlkMiJUJUUUW 
33 33333 J33333J333333 


m 

4> 

O  J  o 


za  at 

UJ«1VJ 

o  o  o 


a  x  x 

UlilAl 
03  O 


A3 
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APPENDIX  D 

INTERPOLATION  FORMULAS 


In  Section  6  we  defined  a  function  B(T)  by  the  infinite  series  (6,6). 
In  order  to  save  computing  tine,  that  function  and  its  derivatives  may  be 
calculated  by  interpolation  in  tables  instead  of  evaluating  the  series*  The 
Interpolation  formulas  are  provided  by  this  appendix. 

the  |uncJiojj  tabljs  in^refer|nc<j»  2,  p.  1144  ff*  contain  the  functions 
B,  Bj,  *  T>  (dB/dT)and  *  T*»(dZB/dT)  •  Hence  in  addition  to  function 
values  at  the  nodes  one  also  has  available  derivatives  of  the  functions. 
This  permits  one  to  use  higher  order  two-point  Hermite  formulas  with 
corresponding  higher  accuracy.  The  basic  Hermite  formulas  are  as  follows. 
Let  the  function  y(x)  and  its  first  derivatives  be  given  at  x  *  0  and  x  *  1, 
and  let  them  be  denoted  by  y0,  yj,  y^*,  and  yj*f  respectively.  Then  an 
approximation  to  y(x)  is 

h(x)  »  yQ  F0(x)  +  y1  Fj(x)  +  y0*  Dq(x)  +  yj*  Dj(x)  .  (D*l) 

The  error  of  the  approximation  is  of  fourth  order  in  x.  The  derivative  of 

(D.l)  is  an  approximation  of  y*(x)  with  a  third  order  error: 


h'(x)  «  y0  F0'(x)  +  yj  Fj'C*)  +  y0'V<x>  +  *1*  Dl’(x)  *  {D*2) 


The  Hermite  factors  In  Eqs.  (D.l)  and  (D.2)  are 


F0(x)  -  (x  -  l)2  (1  +  2x) 
Fj(x)  =  x2  U  -  2  (x  -  1» 
Dq(x)  -  (x  -  l)2  x 
Dj(x)  -  x2  (x  -  1) 


sUttsmaas 


For  an  interpolation  in  the  T-table  between  and  one  has  the  definition 


x  *  (T  -  T  )/(T,-  T  ) 
a  b  a 


and  derivatives  are  transformed  by  using  the  chain  rule 


dx  (Tb  V  d* 


Hence  B  may  be  interpolated  by  the  formula 


*  *  *  *  *  ,Bi-  Biv  . 

B*BaF0  +  BbFl  +  (Tb-Ta)(fiD0  +  *^V  ‘  (D‘6) 

a  b 

*  *  *  * 

In  order  to  coapute  B.  ■  I  1'  one  can  either  Interpolate  B'  making  use 
*  *  * 
of  the  relation  B"  «  B,/t  or  Interpolate  the  function  B  Itself,  or  comp- 

ute  the  derivative  of  Eq.  (D.6),  thereby  disregarding  the  tables.  It  was 
found  from  numerical  experiments  (by  comparing  interpolated  values  with 
exact  values)  that  the  first  approach  is  not  very  accurate.  The  third 
approach  cannot  be  recommended  in  general  and,  therefore,  it  was  not 
tried.  Ihe  final  algorithm  was  based  on  th^  second  approach  and  implemented 
js  follows.  The  derivative  of  the  function  Bj  is  given  the  terms  of  Bj  and 
B^  by  the  formula 


d*  d**  1**  *? 

B.  --(I  B')*jdB'+r 

dT  dT  T 


*  1  *  * 
B">  *  (B.  +  B  > 

T  1  2 


Using  this  value  in  the  interpolation  formula  (D.l)  one  obtains 


*  *  *  ,*  *  ,  ,1a 

B1  -  Bla  F0  +  BlbFl  +  (Tb  -  Ta 


*  *  *  * 

*  *  B,a  +  B-  B  +  B,. 

\  -  V  [  -  —  ■  D0  +  ---  -  V  •  0-8) 

Ta  Tb 


The  function  was  interpolated  by  using  the  inversion  of  Eq.  (D.7): 

*  *  *  *  T  ^®l  * 

Ba’TBi  -Bi-*-^-dT-Bi  *  (D*9) 

b  a 

In  this  formula  the  x-derivatlve  was  coaputed  by  differentiation  of  Eq. 
(D.8),  that  is,  by  Eq.  (D.2).  After  substitution  and  simple  manipulation 
one  obtains 


TbrV 


where  the  relatton  F* 


*  * 

B,  +  B, 
la  2a 


%  *  -  l  *[]  - 

Tb 


has  been  used. 


Limited  comparison  of  the  interpolation^  results  with  exac^t  values 
showed  reasonable  agreement  (eight  digits  for  B,  four  digits  for  B  ).  The 
joutine  UBS  can  easily  be  changed  to  produce  exact  function  values  for  all 
T  by  series  evaluation  instead  of  interpolation.  The  corresponding  increase 
of  computing  time  was  found  to  be  less  than  two  percent. 
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